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ABSTRACT 

The shapelets method for image analysis is based upon the decomposition of localised objects 
into a series of orthogonal components with convenient mathematical properties. We extend 
the "Cartesian shapelet" formalism from earlier work, and construct "polar shapelet" basis 
functions that separate an image into components with explicit rotational symmetries. These 
frequently provide a more compact parameterisation, and can be interpreted in an intuitive 
way. Image manipulation in shapelet space is simplified by the concise expressions for lin- 
ear coordinate transformations; and shape measures (including object photometry, astrometry 
and galaxy morphology estimators) take a naturally elegant form. Particular attention is paid 
to the analysis of astronomical survey images, and we test shapelet techniques upon real data 
from the Hubble Space Telescope. We present a practical method to automatically optimise 
the quality of an arbitrary shapelet decomposition in the presence of observational noise, pix- 
ellisation and a Point-Spread Function. A central component of this procedure is the adaptive 
choice of the shapelet expansion's scale size and truncation order. A complete software pack- 
age to perform shapelet image analysis is made available on the world-wide web. 
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1 INTRODUCTION 

In the shapelets formalism (Refregier 2003; hereafter Shapelets I), 
individual objects in an image are decomposed into weighted sums 
of orthogonal basis functions. The particular set of basis functions 
has been chosen to be mathematically convenient for image ma- 
nipulation and analysis. In astronomical images, it also provides a 
compact representation for the shapes of galaxies of all morpholog- 
ical types. Refregier & Bacon (2003; hereafter Shapelets II) showed 
how these properties could be used to measure the slight distortions 
in galaxy shapes due to weak gravitational lensing. The elegant be- 
haviour of shapelets under Fourier transform also enabled Chang 
& Refregier (2002) to reconstruct images from interferometric ob- 
servations. Massey et al. (2004) used shapelets to simulate realistic 
astronomical images containing galaxies with complex morpholo- 
gies. A classification scheme for galaxy morphologies using Princi- 
pal Component Analysis of the shapelet basis set was discussed in 
that paper and applied to the Sloan Digital Sky Survey by Kelly & 
McKay (2004). A method similar to shapelets has also been inde- 
pendently suggested by Bernstein & larvis (2002; hereafter BJ02). 

In this paper, we expand upon the earlier work of Shapelets I, 
Shapelets II and BI02, developing the formalism of "polar 
shapelets", in which an image is decomposed into components 
with explicit rotational symmetries. Whilst the original Cartesian 
shapelets remain useful in certain situations, the polar shapelets, 
which are separable in r and 6, frequently provide a more elegant 



and intuitive form. We find estimators of an object's flux, position 
and size, that form naturally from groups of its polar shapelet coef- 
ficients. We calculate the behaviour of a polar shapelet model dur- 
ing linear coordinate transformations. We also improve the basic 
shapelet decomposition by incorporating treatments of pixellisa- 
tion, observational noise and point-spread functions, and optimis- 
ing the overall quality of image reconstruction while maximising 
data compression. To test our method upon real data, we extract iso- 
lated galaxies from the Hubble Deep Fields (Williams et al. 1996, 
1998; hereafter HDFs). These deep, high-resolution images from 
the Hubble Space Telescope (HST) provide typical examples of dis- 
tant galaxies' irregular shapes. 

A complete IDL software package to perform the image 
decomposition and shape analyses presented in this paper can 
be downloaded from http://www.astro.caltech.edu/ 
~rjm/ shapelets/. 

This paper is laid out as follows. In §2, we introduce the Carte- 
sian and polar shapelet basis functions, and their relation to each 
other. In §3, we investigate the qualitative effects of varying the 
shapelet scale size /3 and set quantitative goals for the optimisa- 
tion of this choice. In §4, we develop practical techniques to cope 
with the effects of pixellisation, seeing and noise in real data. We 
then demonstrate various applications of shapelets: in §5, we illus- 
trate the manipulation of images in terms of their changing polar 
shapelet coefficients under coordinate transforms. In §6, we con- 
struct basic shape estimators for a shapelet model, including flux. 
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Figure 1. Cartesian shapelet basis functions, parameterized by two integers 
ni and n2, and here truncated at nmax = 6. An image can be decomposed 
into a weighted sum of these functions. This basis is particularly convenient 
for many aspects of image analysis and manipulation commonly used in 
astronomy and other sciences. 



centroid and size measures. In 7, we develop more advanced shape 
measures that can be used to quantitatively distinguish galaxies of 
various morphological types. We conclude in §8. 



2 SHAPELETS FORMALISM 
2.1 Cartesian basis functions 

The shapelet image decomposition method was introduced in 
Shapelets I, and a related method has been independently suggested 
by BJ02. The idea is to express the surface brightness of an object 
f(x, y) as a Unear strni of orthogonal 2D functions, 



oo oo 



m =0 712=0 



"•2 4*fii,n2 /^) 1 



(1) 



where the /m ,n2 are the "shapelet coefficients" to be determined. 
The dimensionful shapelet basis functions 0nin2 are 



<?^ni,n2(x;/3) : 



Hn, (^) (^) 
/3 2"V7rni!n2! 



(2) 



where Hn{x) is a Hermite polynomial of order n, and 13 is the 
shapelet scale size. These Gauss-Hermite polynomials form a com- 
plete and orthonormal basis set; this ensures that the shapelet co- 
efficients for any image can be simply and uniquely determined by 
evaluating the "overlap integral" 



'II 



/(x) 



j(x;/3) d X . 



(3) 



In practice, a shapelet expansion (eq. 1) must be truncated at 
a finite order rii + 112 < nmax. The array of shapelet coefficients 
is sparse for typical galaxy morphologies, which therefore can be 
accurately modelled using only a few coefficients. As shown in 
Shapelets I. data compression ratios as high as 60: 1 can be achieved 
for well resolved HST images. Note however, that our choice of 
Gauss-Hermite basis functions was not governed by the physics of 
galaxy morphology and evolution but by the mathematics of image 
manipulation. As we shall see throughout this paper, a shapelet pa- 
rameterization is mathematically convenient for many tasks com- 
mon in astronomy and other sciences. 

2.2 Polar sliapelet basis functions 

Polar shapelets were introduced in Shapelets I as an orthogonal 
transformation of the Cartesian basis states, and were indepen- 
dently proposed by BJ02. They have all the useful properties of 
Cartesian shapelets, and a similar Gaussian weighting function with 
a given scale size p. However, the polar shapelet basis functions are 
instead separable in r and 6. This renders many operations more 
intuitive, and makes polar shapelet coefficients easy to interpret in 
terms of their explicit rotational symmetries. 

The polar shapelet basis functions 6; ,5) are also pa- 

rameterized by two integers, n and m, and a smooth function 
f{r, 0) in polar coordinates may be decomposed into 



n=0 m=—n 



n,m\n,m 



(4) 



The polar shapelet coefficients /n,m are again given by the "overlap 

integral" 



(5) 



fn,m = I f f{r, 9) Xn,m{r, 9; 13) r drdO . 
J Ju 

BJ02 showed that the "polar Hermite polynomi- 
als" Hni,nr{x), which were described in Shapelets I, are related 
to associated Laguerre polynomials 



Ll(x) = 



pi dxP ^ ' 

for rir > ni by 

H„,,„,(a;) = (-1)"' (n,!)x"'-"'L:^p"'(a:') 



(6) 



(7) 



Here n; and rir are any non-negative integers. In this paper we shall 
instead prefer the simpler n, rn notation, where n = rir + ni and 
m = rir — ni. hi this scheme, n can be any non-negative integer, 
and m can be any integer between —n and n in steps of two. We 
truncate the scries at n < fimax- Although the only allowed states 
are those with n and m both even or both odd, this condition will 
not be written expUcitly alongside every simmiation for the sake of 
brevity. 

As plotted in figure 2, the dimensionful polar shapelet basis 
functions are therefore 



Xn,m{r,0;l3) = 



(-1)- 



2 



(8) 



These are different from the Laguerre expansion used by BJ02 in 
two ways. Those are the complex conjugate of our basis ftmctions: 



Polar Shapelets 3 



-0.4 -0.2 0.0 0.2 0.4 
Polar shapelet basis functions 



2 - 



g 



Real components 



■ 








g 



-6 



" Imaginary components 




1111,111 






















ft 

• 








* 














m 








■ 

ft 



































3 
n 



Figure 2. Polar shapelet basis functions. The real components of the com- 
plex functions are shown in the top panel, and the imaginary components in 
the bottom. The basis functions with m = are wholly real. In a shapelet 
decomposition, all of the basis functions are weighted by a complex num- 
ber, whose magnitude determines the strength of a component and whose 
phase sets its orientation. 



i.e. their m is equivalent to our —m. The Laguerre expansion in 
BJ02 is also normalised by one less factor of /3. This dimensionality 
ensures that, as in the case of Cartesian shapelets, the polar shapelet 
basis functions are orthonormal 



// 

J JM. 



Xn,m{r,0;l3) Xn' ,m'{r,0; P) r drdO = 6„^n'Sm,m' 



(9) 



as well as complete (see e.g. Wiinsche 1998) 

oo n 



n=0 m=—n 



where 5 is the Kronecker delta and the asterisk denotes complex 
conjugation. Only those basis fimctions with m = contain net 
flux. 



(11) 



/ / Xn,m(r, e- P) r ArA6 = 2^(5 5^o ■ 
J Jr 

Figure 3 demonstrates the polar shapelet decomposition of a 
galaxy found in the HDF. The original image (middle left panel) 
agrees well with the reconstruction using rimax = 20. The top 
panel shows the modulus of the polar shapelet coefficients as func- 
tion of the n and in indices. The dominant coefficients have small 
values of both indices, demonstrating the compactness of a polar 
shapelet representation, and further improved prospects for data 
compression. The bottom panel shows the reconstruction of the 
galaxy using only coefficients with given values of n or \m\, thus 
highlighting the contributions of terms with specific rotational sym- 
metries. The off-central bulge is captured in the \m\ = 1 coeffi- 
cients, and the main spiral arms in the \m\ = 2 coefficients. The 
spiral arms can also be seen as the rotation of the n-only recon- 
structions with increasing radius. The fainter spiral arms appear as 
an interplay of the |m| =4 and 5 coefficients. 



2.3 Conversion between Cartesian and polar shapelets 

Cartesian shapelets are real functions, but polar shapelet basis func- 
tions Xn,m. and coefficients fn,mWt both complex. However, their 
symmetries 



Xn,-m(r, 6; 0) = xl.m{r, 9; 0) = Xr^,m{r, -6; f3) , 



(12) 



simplify matters somewhat if we are concerned only with the rep- 
resentation of real functions /(x), like the surface brightness of an 
image. Equations (9) and (12) imply that /(x) is real if and only if 

fn. — va = fn,m • (13) 

Coefficients with m = are thus wholly real. All polar shapelet 
coefficients are paired with their complex conjugate on the other 
side of the line m — 0. Therefore, even though the polar shapelet 
coefficients /„ m are generally complex, the number of indepen- 
dent parameters in the shapelet decomposition of a real function is 
conserved from the Cartesian case. 

A set of Cartesian shapelet coefficients fm ,n2 with ni+n2 < 
^max can be transformed, into polar shapelet representation with 
n < rimax, using 

f _ n-^^m ni\n2\ ^ 



n-'r'm \ | / n — m \ \ 
K 1 }-\ 2 )' 



E E 



V 2 ] 



n — m \ 
, 2 ) 

n'l 



The particular choices of truncation scheme for Cartesian and polar 
shapelets now make sense as a way to keep this mapping one-to- 
one. 



3 CHOICE OF SHAPELET SCALE SIZE 

A shapelet decomposition requires values for the scale size /3 and 
for the centre of the basis functions Xc to be specified in advance. 
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Figure 3. Example polar shapelet decomposition of a HDF galaxy. Top 
panel: the moduli of the polar shapelet coefficients, with a logarithmic 
colour scale. Bottom panel: the original galaxy image using a linear colour 
scale and its shapelet reconstruction using Umax = 20. Additional recon- 
structions are shown using only particular sets of coefficients, to highlight 
the contribution of components containing different symmetries. 



Choosing the centre is relatively easy: there are many methods 
well-known in the astronomical literature to accurately determine 
astrometry from the flux-weighted moments of objects. However, 
the selection of /3 is a less well-posed problem. In this section, we 
shall first use some properties of polar shapelets to describe the 
effect that the choice of the scale size has upon a shapelet decom- 
position. We shall then set quantitative criteria for the selection of 
P in arbitrary galaxy images that we can implement in a practical 
algorithm. 

Note that the selection of /3 will be linked to the selection of 
Wmax- As shown in Shapelets 1 §2.4, these two parameters deter- 
mine the maximum extent 6'inax and finest resolution ^min that can 
be successfully captured by a shapelet model. If Umax — > oo, any 
object can be represented using any scale size B. But if the shapelet 
expansion is truncated at finite rimax, the shape information needs 
to be more efficiently contained within fewer coefficients. It is 
clearly desirable in this situation to select a scale size /? that corn- 



Table 1. The first few rotationally-invariant polar Shapelet basis functions. 



X6.oM=^[l-3^ + |^-i^]e^ 
X8,o(r)= ^ [l - 4^^ + f ^ - 1^ + 

XlO,o(r-)= ^[l-5f^ + f^-f§^ + Af^ + Tiof^]e^ 

presses information, and lets us store the smallest possible number 
of coefficients. 

3.1 Radial profiles 

Our discussion can be simplified by initially considering only 
the radial profile of an object, thus reducing the task to a one- 
dimensional problem. Let us consider an object with surface bright- 
ness /(x). The object's radial profile f(r) is its brightness averaged 
in concentric rings about its centre, i.e. 

l{r) = ^ fj f{r,e) AO. (15) 

With the object decomposed into polar shapelets as in equation (5), 
it is easy to show that this is given by 

even 

71 

This simple expression results from the fact that only the m = 
basis functions are invariant under rotations. These are given by 

X„o(r;/3) = tir^L%[r' j (3')e-^' . (17) 

The first few rotationally-invariant basis functions are written ex- 
plicitly in table 1. 

As a concrete example, we consider the decomposition of 
galaxy images from the Hubble Deep Fields (Williams et al. 1996, 
1998). The mean radial profile of spiral galaxies is typically expo- 
nential, /(r) oc e'"^^^" , with some characteristic radius scale ro. 
Figure 4 shows the shapelet reconstruction of an exponential radial 
profile using various values of /3, with rimax = 20 and the integral 
in equation (5) calculated numerically. 

As can be seen in the top panel of figure 4, the quality of 
the reconstruction depends on the choice of fl. For small values 
{f5 ^ 0.4ro) the reconstruction is oscillatory and cuts off the profile 
at large radii (r > 1.5ro). For large values (/? ^ 0.8ro), the recon- 
struction fails to reproduce the cusp at small radii (r ^ 0.4ro) and 
exceeds the true profile at r ~ 0.6ro. However, for intermediate 
values (O.Sro < < l.lro), the reconstruction is good through- 
out the range O.lro ^ r ^ 2.8ro. This range can of course be 
expanded by including more shapelet coefficients of higher order. 
As Timax —> oo, the input model can be recovered with arbitrary 
precision using any value of 0. 

The corresponding behaviour in shapelet space is apparent in 
the bottom panel of figure 4. The fn,o coefficients can be thought as 
the profile of the galaxy in shapelet space or the "shapelet profile". 
For low values of /? the shapelet profile is very flat, showing that 



Polar Shapelets 5 





10 

n 



15 



20 



Figure 4. Decomposition of an exponential profile into radial polar 
shapelets. Top panel: The thick dark line shows the input exponential pro- 
file. The reconstructed profile is shown for different values of the shapelet 
scale j3 with Wmax = 20. Bottom panel: the corresponding shapelet coeffi- 
cient profile fnQ versus shapelet order n. 



the power is distributed almost evenly throughout all orders. For 



(5 = 0.5ro, the coefficients 



are seen numerically to be oc 



(n + 1)" . This will be an important result for the convergence of 
shape estimators formed from series of shapelet coefficients in §6. 
Convergence is fastest at /3 « 0.8ro, with On^ni oc (n + 
For higher values of /3, the signs of an.m begin to alternate, and the 
convergence once more falls below oc (n + at /3 « l.lro. 

Figure 5 demonstrates the importance of a proper choice of 
the parameters /3 and rimax for the practical decomposition of spi- 
ral galaxy in the HDF. Its spiral arms complicate measurement, but 
its radial profile is approximately an exponential with a scale length 
of ro ~ 3" (12 pixels). The left column shows the growth in com- 
plexity of a shapelet model using increasing rimax- Note in particu- 
lar the rotation of the core ellipticity as rimax is increased from 2 to 
8 and higher order moments are included to resolve the spiral arms. 
In this column, 13 is allowed to vary in order to minimise the least- 
squares difference between the model and the HDF image, shown 
at the bottom. The middle column shows shapelet decompositions 
at fixed rimax = 20, with varying 13. The residuals are plotted in 
the right hand column. As in figure 4, we find that the best overall 
image reconstruction uses 0.5ro ^ /3 ^ 0.7ro. This is perhaps at 
the low end of the range suggested by figure 4 because of the extra 
high-frequency detail contained in the spiral arms. 

By experimentation we have found a fairly wide range of (3 
values that produce a faithful shapelet reconstruction. The informa- 
tion is then concentrated into the few lowest shapelet states, with 
fast convergence to the final model, and truncation is possible at a 
computationally acceptable value of rimax- We shall now consider 
ways to formalize this process, and hone our choice of Xc, P and 
rimax using quantitative criteria. 
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Figure 5. Shapelet decomposition of a real spiral galaxy in the HDF. The 
best-fit de Vaucouleurs profile has rg ~ 12 pixels. Left-hand column: the 
shapelet model shows growing complexity with increasing rimax. For each 
of these fits, (3 is varied to minimise the least squares difference between the 
data and the model. Right hand columns: the shapelet decomposition has a 
preferred scale size. The residual between the original (in the bottom-right 
panel) and these models with fixed rimax = 20 and varying /3, is smallest 
with (3 ~ 0.5ro. 
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3.2 Existing optimisation methods 

Methods in the literature that face similar choices suggest several 
distinct philosophies for the quantitative selection of parameters 
equivalent to Xc, P and rimax- The suggestions, outiined below, dif- 
fer both in the goals set for for an ideal decomposition and the 
method used to achieve it. 

• Shapelets I suggested a geometrical argument using 6'jnin, 
^max: the minimum (PSF or pixel) and maximum (entire image) 
sizes on which information exists. This could be iterated using 
functional rules on Xc and r] as defined by shapelet coefficients. 
However, the coefficients change as a function of nmax, and it is 
not clear what the rules should be. 

• Van der Marel & Franx (1993) fit ID Gauss-Hermite poly- 
nomials to spectral lines. They arbitrarily fix rimax = 6, prob- 
ably finding this sufficient because their spectra have relatively 
high S/N and their lines have a nearly Gaussian profile. Parameters 
equivalent to Xc and (5 are obtained from the best-fitting Gaussian. 
This also determines /o and in ID is equivalent to constraining 
h = h = 0, i.e. the derivatives of the Gaussian with respect to 
Xc and p. The number of variables is reduced and the problem ren- 
dered tractable. Unfortunately, this does not help us in 2D because 
while both ai±i can be forced to zero by varying Xc, no unique 
recipe can be found for setting the three n = 2 states using only 
one value (3. 

• Van der Marel et al. (1994a,b) relax the constraint on fi. This 
is an improvement since /i is only the first term of an expression 
for the centroid, expanded using all odd /„ in equation (51). With- 
out the higher order corrections, the true object centroid is moved 
slightly from the origin: amongst other things rendering rotations 
and shear operations more complicated. Instead, they set Xc from 
the theoretical rest wavelength of a line. Unfortunately, astrometric 
calibration clearly cannot be done with such accuracy. Nor has the 
n = 2 problem been solved. 

• Kaiser, Squires & Broadhurst (1995) combine fitting with a 
stand-alone object detection algorithm, hf indpeaks. Translated 
into shapelet language, their approach is roughly equivalent to plac- 
ing Xc at data peaks then finding a width /3 such that signal-to-noise 
ratio V in /oo is maximised. 

• Bernstein & Jarvis (2002) propose a similar approach. They 
prescribe /3 by requiring /20 = 0, while moving x,, to ensure 
/i±i = 0. Higher coefficients are then determined afterwards by 
linear decomposition. To first order, this (3 constraint is equivalent 
to that for hf indpeaks. This /3 is generally larger than values 
chosen by our method below, and it can be several times larger 
for a high signal-to-noise object containing lots of substructure like 
the galaxy in figure 3. This method may indeed prescribe the op- 
timal decomposition for weak lensing as the shear signal in the 
quadrupole moments becomes concentrated in one number; how- 
ever, a predisposition towards particular states often leads to poor 
overall image reconstruction and PSF deconvolution, so it is not 
necessarily ideal for all applications. 

• Kelly & McKay (2004) were able to set a fixed physical scale 
of =2kpc for galaxies in the Sloan Digital Sky Survey, where 
photometric redshifts were available. However, galaxies have a 
broad distribution of physical sizes, and it may in fact become more 
difficult to interpret a shapelet model derived using this method. 

• Marshall (in preparation) describes a fully Bayesian approach 
to applying the shapelet transform in the context of image recon- 
struction. Here, Xc, (3 and rimax are varied in order to maximise the 
evidence (the probability of observing the data, marginalised over 
all shapelet coefficients). At high S/N, this method gives a value of 



(3 which approaches the same as that from our method below, 
but otherwise tends to prefer a fractionally larger (3, conservatively 
eliminating some 'noise' in favour of a smoother image reconstruc- 
tion. However, this is computationally slow, a serious issue when 
analysing large images. 

3.3 Optimisation of image reconstruction 

We shall adopt a choice of and rimax that is suitable for many 
applications, including overall image reconstruction and PSF de- 
convolution. Different models will be quantitatively compared via 
the overall reconstruction residual 

2 _ (/obs(x) - /,ec(x;/3))^\/-i(/obs(x) - /rec(x;/3)) 

Xr = — I (to) 

ripixels ricoefFs 

where /obs(x) is the observed image, and /rcc(x; i'?) the recon- 
structed image from the shapelet model. V is the covariance matrix 
between pixel values, i.e. its diagonal elements are the noise vari- 
ance in each pixel. We will need to know this a priori, or estimate it 
from blank areas of the image, ripixeis is the number of pixels in the 
observed image, and rrcoofts is the number of shapelet coefficients 
used in the model. The residual itself has variance (Lupton 1993) 

'y'\xi) = z ■ (19) 

^pixels ricoeffs 

An example of typical isocontours on an rimax (3 plane 
is shown in figure 6 for an elliptical galaxy from the HDF. The 
horizontal trough is present for all galaxies (and many other iso- 
lated objects). This demonstrates that there is indeed an optimum 
(3 for the reconstruction of this image. As one might expect, it is 
roughly independent of rimax, but decreases very slightly as more 
coefficients are added. By increasing rimax — > cx3, the reconstruc- 
tion can be improved to arbitrary precision. However, stopping at 
the Xr = 1 contour produces a model whose residual is consis- 
tent with the noise. Additional coefficients would just model the 
background noise and should be excluded. 

The form of these typical contours thus suggests a unique lo- 
cation in parameter space. We will choose (3 and rimax so that the 
model lies at the intersection of the trough and the — \ contour, 
i.e. at the left-most point on the contour. To achieve this, we set 
quantitative goals of 

dp 



0, 



Xc =0. 



(20) 
(21) 



and 



Xr 



1 or flattens out 



dxl 



< <y {xl) 



(22) 



The first constraint ensures that the scale size is well-suited to ef- 
ficiently model the image. The second ensures that the shapelet 
center matches the object centroid. The third guarantees that suffi- 
cient coefficients (rimax) are included to model an object, but with 
truncation that 'smooths over' observational noise. A flatness con- 
straint is also included (in the right hand side of equation 22). This 
is particularly important for galaxies with a near neighbour or for 
very faint objects that have noisy and fragmented Xr contours. In 
these cases, including additional shapelet coefficients may not sig- 
nificantly improve a fit, so the series is truncated early. 

We apply extra geometrical constraints to the minimum ^min 
and maximum 0max scales of the decomposition, to prevent the 
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Figure 6. Xr isocontours on an rimax 'os (3 plane for an elliptical galaxy in 
the HDF. The roughly horizontal trough is typical, with a well-pronounced 
excursion of the Xr contours to lower Wmax for well-chosen values of /3. 
The challenge is to locate the leftmost section of the Xr = 1 contour in 
an automated and efficient way. The arrows show individual steps (each 
containing several sub-steps) taken by our optimization algorithm described 
in §4. Also shown are geometrical Snjinj^max constraints and the target 
X? = 1 contour. 



model from containing features smaller than the pixel scale or ex- 
tending off the edge of an image, where it would be unconstrained. 



3.4 Automatic optimisation algorithm 

Satisfying the three conditions (20), (21) and (22) would ensure 
that a shapelet decomposition uses the optimum values of rimax, 
/3 and Xc. It is easy to determine the values of these parameters 
once the entire rimax vs /3 plane has been examined, as in figure 6. 
However, this is a slow process, so we need a practical algorithm to 
more efficiently explore this parameter space, and to iterate rapidly 
towards our targets. The numerical implementation of this iteration 
will inevitably be non-trivial, because it combines both minimisa- 
tion and root finding, in a space with one axis discrete. Here we 
describe a code that we have developed to repeatedly decompose 
an object into shapelets, test the residual, and improve the decom- 
position parameters. Its stepwise approach is shown in figure 6, and 
the full code can be downloaded from the world wide web. 

Objects are first detected in an image using SExtractor 
(Bertin & Arnouts 1996), a friends-of-friends peak-finding algo- 
rithm. After experimenting on various data sets, we have found the 
results of SExtractor highly sensitive to input settings. To avoid 
reliance upon these setting, we use SExtractor as sparingly as 
possible. We set low detection thresholds in order to obtain a com- 
plete catalogue, and filter out false detections later. We use the mea- 
surement of each object's FWHM to make an initial guess at /3, and 
also to set the size of the fixed, circular "postage stamp" region that 
is extracted around each object. We aim for a postage stamp large 
enough to contain the entire object, but small enough to isolate it 
from its neighbours and to make the routine computationally effi- 



cient. We then use the SExtractor segmentation map to identify 
pixels in the postage stamp but well away from any object. These 
are used to estimate the background noise level, or to locally renor- 
malise the pixel weight map. Within reasonable limits, the process 
is stable with respect to such parameters and we shall not be too 
concerned as to the exact SExtractor settings. 

Using constant Wmax = 2 for speed, 13 is varied in order to 
minimise Xr and satisfy the criterion in equation (20), via a ID 
version of the Numerical Recipes AMOEBA routine: crawling verti- 
cally in figure 6. During each step of this iteration, the centroid is 
simultaneously shifted to re-zero the series in equation (51) in the 
shapelet coefficients and thus satisfy the criterion in equation (21). 
Since the calculation of the centroid is independent of P for iso- 
lated objects (see §7), this part of the iteration is both stable and 
fast. Figure 6 also shows the additional geometrical constraints of 
6min > 0.2 pixel and 6'nmx not falling off the edge of the postage 
stamp. These act as hard boundaries to the region of parameter 
space that the amoeba is allowed to explore. 

Once the optimum (3 has been found, nmax is increased until 
the criterion in equation (22) is satisfied: crawling horizontally in 
figure 6. The increases are done in steps of two, because even n 
states frequently improve the fit more than odd n states (primarily 
due to the additiona of a new /„,o circular state). The value of nmax 
is fine-tuned to the exact best value at the end. If two values of nmax 
both allow a decomposition with y^. — \ ± la, the lower value is 
taken. 

If the object warranted more coefficients than the initial guess 
of nmax = 2, /3 and Xc are again readjusted at the new rimax, using 
our ID AMOEBA routine. Another nmax search then starts back 
at nmax ~ 2 and increases again in steps of two. The algorithm 
terminates when either the horizontal or vertical search returns to 
the value it started with. All three conditions in equations (20), (21) 
and (22) have then been met. Computation time for each object 
increases oc nmax - On a single 2Ghz processor, our algorithm takes 
about 45 minutes to process all of the 3596 objects detected in the 
HDF North. 

A selection of reasonably bright HDF galaxies is shown with 
their shapelet models in figure 7. The right-hand column shows 
the reconstruction residuals, which are consistent with noise even 
for irregular galaxy morphologies. A comparison of their shapelet- 
based shape estimators to traditional SExtractor measurements 
is shown in figure 8. 



4 DECOMPOSITION OF REAL DATA 
4.1 Least squares fitting 

Unlike the continuous, analytic formalism presented in §2, real im- 
ages are complicated by pixellisation, PSF convolution and noise. 
In order to incorporate these effects, we shall first adopt a some- 
what different approach to shapelet decomposition than the overlap 
integrals (3) and (5). We shall instead fit shapelet coefficients to the 
data using a least-squares method. Since the model /rcc (x) in equa- 
tion (18) is linear in the shapelet coefficients, we can solve for the 
minimum Xr solution (18) exactly. We obtain (see Lupton 1993; 
Chang & Refregier 2002) 



(23) 



where fn,m is a vector of the derived shapelet coefficients, fx,y 
the surface brightness in each pixel arranged as a data vector, V 
the covariance matrix between pixel values and M is a matrix of 
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Real HDF image Shapelet model Residual 



i 



Figure 7. Shapelet models of a selection of HDF galaxies, with their 
shapelet scale size /3 and maximum order nmax determined automatically. 
In all cases, the image residuals are entirely consistent with noise. Our code 
to perform this task, by minimising the least squares difference between the 
model and input images, is described in the text. 
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Figure 8. The successful recovery of object statistics from the shapelet pa- 
rameters of HDF galaxies. For comparison to the SExtractor measure- 
ments, shapelet size measurements are shown without PSF deconvolution. 
In the right-hand panels, galaxies requiring nmax > 15 coefficients have 
been forced into the final bin, and in the bottom-right panel, points have 
been randomly offset a small amoimt for clarity. 



each shapelet basis function evaluated in each pixel. A fit achieving 
Xr = 1 has successfully modelled all significant spatial variation 
in the image, and removed observational noise. 

If the noise per pixel is known, Icr confidence limits can be 
derived on all of the assigned coefficients using this fitting method 
(Lupton 1993). If a complete pixel noise map is available {e.g. from 
multiple exposures stacked using DRI Z ZLE software - Fruchter & 
Hook, 2002), it can be used to down-weight noisy pixels where 
cosmic rays or hot/cold pixels were present in some of the expo- 
sures. Although the code available on the world-wide web simply 
uses a diagonal matrix for V that contains only the noise level in 
each pixel, the method is, in general, able to use the full covariance 
matrix that contains the amount of covariance between different 
pixels. In real data, the flux in adjacent pixels is indeed slightly 
correlated because of convolution with the PSF and also because of 
additional aUasing effects introduced by DRIZZLE. If this effect is 
important, the pixel-to-pixel covariances could be estimated from 
empty regions of an image and included in the calculation. In par- 
ticular, this may have a small improvement on statistics measured 
from very small objects (c./ Massey et at. 2004). 

A constant background level can also be removed using this 
method, by adding an undetermined constant to the set of basis 
functions. Poor flat fielding or local background gradients near a 
bright object can also be fit and removed by adding a plane with 
variable slope. Although these functions are not strictly orthogonal, 
the procedure works well in practice as long as there are sufficient 
pixels around the fitted object that contain only background noise. 



4.2 PSF deconvolution 

All real images are inevitably seen after convolution with a Point 
Spread Function (PSF). In astronomy, this is typically caused by 
atmospheric turbulence or "seeing" (for ground-based observato- 
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ries), aperture diffraction at tlie primary mirror, and imperfect tele- 
scope tracking or optics. The combination of such effects can be 
measured from the size and shape of stars observed in an image 
(because these distant objects would be point-like in the absence of 
a PSF), and can be fit with a shapelet model in the same way as the 
galaxies. Shapelets I presented the matrix operation for convolv- 
ing an image with a Gaussian PSF in shapelet space. Shapelets II 
extended this derivation to a general PSF and demonstrated PSF 
c/econvolution via matrix inversion. However, the inversion of the 
PSF matrix is potentially slow and may be numerically unstable. 
Our least-squares fitting method will allow us to elegantly sidestep 
this process by convolving the basis functions with the PSF model 
in advance, then fitting this new basis set to the data. The returned 
shapelet model, reconstructed using the unconvolved basis func- 
tions, will be automatically deconvolved from the PSF. 

The formalism for convolution in shapelet space is presented 
in Shapelets I §4 and involves three separate scale sizes for three 
separate objects: a for the unconvolved model, for the PSF, 
and 7 for the convolved model (there are also corresponding val- 
ues of nZ. 
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). We assume that ,3 is known. We Ti^Tim plus Shapelet model, 



can optimise a as in section §3.3. However, the choice of 7 is a 
matter entirely internal to the fitting procedure. Just as before, if 

J^max 00, any 7 will work (but this time without increasing the 
number of external free parameters in the model). In practice, how- 
ever, it is still necessary to truncate this series somewhere. Note 
that 7^ = _^ (3^ was incorrectly suggested as a "natural choice" 
for this parameter in Shapelets I. Another choice would be 7 = a, 



Residual 



which, with 



makes the convolution matrix Pn. 



symmetric and thus simplifies its calculation. 

The optimum values for 7 and nJiax aie in fact obtained from 
an argument concerning the information present in shapelet coeffi- 
cients. A shapelet model contains information only between a min- 
imum and maximimi scales 



+ 1 



and On 



During convolution, O'j^i^ and O'^^^ add in quadrature to produce 



. + 1 



(24) 



e^,^ and eL 



'mm' -max — -max add similarly to produce OZie^; the range 
of scales on which information is available decreases or remains 
constant. Values of 7 and n^iax can be chosen to most efficiently 
capture the information contained between the new scales. Writing 

(^max + 1) as A^Q, etc. for brevity, we find 



7 = 
and 



''max 



a^Ng + B^Nc 



NcNp - 1 



(25) 



(26) 



The PSFs of cameras on board the Hubble Space Tele- 
scope are well known and stable. Figure 9 shows an oversampled 
TinyTim (Krist 1997) model of the Wide-Field Planetary Cam- 
era 2 (WFPC2) PSF, raytraced through an engineering model, plus 
charge diffusion to simulate photon capture within the CCD cam- 
eras. This is easy to model with shapelets, except for the fact that 
its steep cusp and extended wings are intrinsically ill-matched to 
the Gaussian around which shapelets are constructed. The PSF is 
shown in the figure beside a shapelet decomposition up to rimix = 
15. This is sufficient to accurately capture the core and the first 
two diffraction rings, which are already more than two orders of 
magnitude below the maximum, but does not extend to the four 
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15 



Figure 9. Shapelet model of the TinyTim (Krist 1995) WFPC2 PSF plus 
charge diffusion. Top panel: a horizontal slice through the centre of the 
PSF. Bottom panel: the moduh of its polar shapelet coefficients to nmax = 
15. Note that the amplitude scales are all logarithmic: the core is actually 
modelled very successfully out to the second diffraction ring. For speed we 
do not bother capturing the wings. 



faint diffraction spikes or far into the low-level wings (note that 
the colour scales are logarithmic). In principle, this could be fur- 
ther extended at a cost to processor time by using more shapelet 
coefficients. 

Figure 10 demonstrates successful PSF deconvolution. A 
galaxy from the HDF is convolved with the WFPC2 PSF (in real 
space). This is treated as the observed image, and deconvolved from 
the PSF using a shapelet fit. The resulting reconstruction is in good 
agreement with the original galaxy image, as can be seen from its 
small residual. Note that the optimum scale size /3 for the model is 
slightly lower when PSF deconvolution is performed. This reflects 
the need to capture finer details. 

4.3 Pixellisation 

Real image data is typically stored in discrete pixels. To Unk this to 
the analytic shapelet formalism, one must either smooth the data or 
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Figure 10. Demonstration of deconvolution from an observational PSF. 

Top-left panel: a real HDF galaxy. Bottom-left panel: the WFPC2 PSF 
model, from figure 9 but displayed here with a hnear colour scale. Bottom- 
middle panel: the galaxy image convolved with the PSF. This convolution 
has been performed in real space, to disassociate the operation from any- 
thing involving shapelets. Top-middle panel: a shapelet reconstruction and 
deconvolution of the galaxy to rimax = 20, obtained from a fit to the con- 
volved image, assuming knowledge of the PSF. Top-right panel, the differ- 
ence between the true galaxy image and the shapelet model after deconvo- 
lution. This small residual demonstrates the success of PSF deconvolution 
using shapelets. 



pixellate the shapelet basis functions. Smoothing the data requires 
an arbitrary interpolation scheme to be defined, and resampling the 
data onto smaller pixels can be very slow. A better approach is to 
leave the data alone, and discretize the smooth shapelet basis func- 
tions. This reduces the integrals in equations (3) and (5) to sums 
over pixel values, which are fast to compute. However, they are no 
longer analytically exact. We therefore need to define a discretiza- 
tion scheme that keeps the basis functions as orthogonal as possi- 
ble, and the integrals as accurate as possible. 

As pointed out by Berry, Hobson & Withington (2004), one 
cannot simply adopt the value of basis fimctions at the centre 
of each pixel. Basis fimctions that contain oscillations on scales 
smaller than the pixel size are sampled in an essentially random 
manner. Their discrete versions are then neither representative of 
the analytic function nor orthogonal. Degeneracies are introduced 
between shapelet coefficients during the decomposition that un- 
stabilise the inversion of coefficient matrices in the reconstructed 
model, and bias quantities like an object's flux. Fortunately, this is 
rarely a problem in practical cases, because we can choose rimax 
and /3 in advance to isolate only those basis functions that contain 
oscillations on scales larger than the pixel (or seeing) size. Under 
these conditions. Berry, Hobson & Withington (2004) show that the 
shapelet basis functions are indeed orthogonal. 

We suggest an even safer alternative here. The Cartesian basis 
functions are separable in x and y, and may be analytically inte- 
grated within rectangular pixels. This is exactly the same process 
undergone by photons arriving at a CCD, where the smooth func- 
tion of a real scene gets binned into digital squares. Once we have 
convolved the basis functions with the PSF, and integrated them 
within pixels, they can be suitably matched to the data. 



To integrate the 2D Cartesian basis functions, first consider the 
ID basis functions from Shapelets I, 



(27) 



Integrating by parts and using two well-known identities (see e.g. 
Boas Ch. 12) 



Hn{x) = 2xHn-i{x) - 2{n - l)i?„-2(a;) 
and 

dHn-i{x) 



dx 



2(n- l)i^„-2(x) , 



one can obtain the recurrence relation 



dx 



In = I (j>n{x) 

/2 r lb 
-Un-l(x) + 
n L J a 

Finally, note that 



n- 1 



lo = 



Pn2 
2 



erf 



and 



(28) 
(29) 

(30) 
(31) 

(32) 
(33) 



This supplies all the necessary integrals. Since the 2D Cartesian 
basis functions are separable in x and y, it is easy to extend this 
derivation to integrate within square CCD pixels: 



J J ' 



L {x)4>n2 iv) Ay = Ini X 7„ 



(34) 



where, if there is no 'dead zone' aroimd the edge of a pixel, (61 — 
oi ) X (62 — 02) is the angular size of a pixel. A missing pixel border, 
due for instance to electronics which is unresponsive to light, can 
be included by altering the limits on the integral. 

We can either use this result to obtain a model in Cartesian 
shapelet space, which can later be converted to a polar shapelet 
representation using equation (14), or we can integrate the polar 
shapelet basis functions within pixels using the same equations. 
This integration is a particularly important advance for small galax- 
ies or for shapelet basis functions at high-n, that can contain oscil- 
lations smaller than a single pixel. 

The symmetries of polar shapelets can also be used to integrate 
models within circular apertures using equations (46) to (49). 



5 COORDINATE TRANSFORMATION 

Image manipulation via linear transformations is simple in shapelet 
space. As in Shapelets I, let us consider an infinitesimal coordinate 

transformation x ^ (1 + ^)x + e, where e = {ei, €2} is a dis- 
placement and ^ is a 2 X 2 matrix parametrized as 



K + 71 
72 + P 



72 
K - 



- P 
71 



(35) 



The parameters p, k, e and 7, correspond to infinitesimal rotations, 
dilations, translations and shears. 

An image transforms as /(x) — » /'(x) ~ /(x — *x — e), 
which can be written as 



/' ~ (1 -h + /t^ + ijSj -h eifi)f, 



(36) 
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where R, K, Si and Ti are the operators generating rotation, con- 
vergence, shears and translations, respectively. We adopt a notation 
from weak gravitational lensing, where a "convergence" k corre- 
sponds to a change in an object's radius by a factor (1 — k)~^. 
These transformations can be viewed as a mapping of .fn,m coeffi- 
cients in shapelet space. For example, a finite rotation is 



,71 



J n, 



SO a rotation through 180° can be written as 

-Rl80° : fn,m f'n^m = (" 1)™ fn,m ■ 



(37) 



(38) 



An (infinitesimal) dilation can be performed in polar shapelet 
space by mapping the shapelet coefficients as 



K:fn 



f 



(1 + «) /n 



(39) 



+ 



2 \/ {n - rn) {n + m) fn-2.-n 



jn + 2){n + jn + 2) fn+2,, 



The shapelet model may require more coefficients after this trans- 
formation. Note that this dilation operation increases both the flux 
and the image area by a factor 1 + 2k, thus conserving surface 
brightness. To instead perform a dilation that conserves the total 
flux, divide the right hand side of equation (39) by this factor. To 
first order, this is 



K : fn,m *■ fn,m — (1 fn,; 



(40) 



+ '^^/{n-m){n + m) /n-2,m 

- ^\/{n- m + 2)(n + m + 2) f„+2,m. ■ 

In §6, we shall ensure that shape estimators for a shapelet model 
are independent of the scale factor chosen for the decomposition 
by ensuring that the estimators are unchanged under this mapping. 

Rather than these first-order approximations, finite dilations 
can be performed to all orders using the rescahng matrix in the 
appendix of Shapelets I. This is identical to the convolution matrix, 
but the image is convolved with a (5-function. 

Shears and translations can be performed using 



S '. fn,m ^ fn,m — fn, 
71 + *72 



(41) 



+ - 



and 

f:f„ 



I V^(n + m)(n + m - 2) /„-2,m-2 

- ^/{n-m + 2)(n - m + 4) /„+2,m-2| 
- '^^ ^'^'^ I ^J{n- ■m){n - m - 2) /„_2,m+2 

- \/(n -h m -I- 2)(n -h m -I- 4) /„+2,m+2 j 

*■ fn,m — fn,Tn (42) 
+ ^^^^ {V{n + m) fn-l,m-l 

- ^/{n-m + 2) /„+i,m_i j 

I l/(n - m) /n-l,m+l 
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with the translation specified in units of (3. 
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Figure 11. Some simple operations applied to a real galaxy image, by using 
the polar shapelet ladder operators or coefficient mappings as described in 
the text. The central image is the original galaxy. Starting at the bottom-left 
and proceeding clockwise, the other images show rotation by 40°, dilation 
of K = 0. 15, shears of 7 = 10%, translations, circularisation and reflection 
in the a;-axis. 



Other image manipulations can also be represented as map- 
pings of shapelet coefficients. Changes of flux by a factor B are 
trivially implemented by the mapping 



B : fn,m * /n, 



B X fn,, 



(43) 



It is also possible to circularise an object with the mapping (see 
§3.1) 

C '. fn,m *■ fn,m ~ fn,m ^mO ? (44) 

or to flip an object's parity by reflection in the a;-axis using 

P ■ /n,m *■ fn,m — fn,m ■ (45) 

Combining this P with the rotation operator allows reflections to 

be performed in any axis. 

These operators' actions are demonstrated upon a real galaxy 
image in figure 11. 



6 OBJECT SHAPE MEASUREMENT 

The above symmetries of the polar shapelet basis functions can be 
used to identify combinations of shapelet coefficients that measure 
an object's flux (photometry), centroid position (astrometry) and 
size. Similar weighted combinations of Cartesian shapelet coeffi- 
cients were found in Shapelets I, but we find the interpretation of 
polar shapelets more intuitive, and the expressions below are usu- 
ally more simple than their Cartesian equivalents. For example, the 
rotationally invariant part of an object is isolated into its m = 
coefficients. The linear offset of an object from the origin is de- 
scribed by its m = ±1 coefficients and the ellipticity of an object 
by its m = ±2 coefficients. In the latter cases, the magnitude of 
the coefficients indicate an ampUtude, and the phases a direction. 
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6.1 Photometry 

Practical measurements of objects' flux usually introduce a Gaus- 
sian or top-hat weight function in order to limit contamination from 
surrounding noise and nearby objects. The flux of a shapelet model 
inside a circular aperture can be calculated using only the coeffi- 
cients with m = 0. All other coefficients correspond to basis func- 
tions with positive and negative regions that cancel out under in- 
tegration around 6. From equations (16) and (17) for an object's 
radial profile, we find that 



6.2 Astrometry 

It can similarly be shown that the unweighted centroid {Xc, Vc) is 
/,4(x + ij/)/(x)d^x ifint^p^"^" 



Xc+iyc = 



^(n+l)i/„i(51) 



Here the summation is over only odd values of n, because only 
these have the m — ±1 coefficients that possess the desired rota- 
tional symmetries. 



,.2n i-R ^ even 

/ / /(r) r drd^ = (47r) ^ /3 V /„.o In 
Jo Jo 



(46) 



where 



g 2^ r dr . 



(47) 



Using relation (A9) to integrate by this parts, we can find a recur- 
sion relation* 



In = (-1)5 i 1 - 4 ) - 2 E \ , (48) 



and a closed form 



/„ = 1 - 



i/2 



(49) 



However, the imposition of a integration boundary is unnec- 
essary with shapelets because the model is analytic and noise-free. 
In the limit of R ^ oo, we obtain a simple expression for the total 
flux in a shapelet model 



= // /(x) <fx = (47r)i/3 fno , 
J Jm 



(50) 



a result that can also be recovered by transforming the sum over 
Cartesian shapelet coefficients from Shapelets I into polar shapelet 
space via equation (14). Cartesian shapelet models can also be in- 
tegrated within square apertures using equations (30) to (34). 

This extrapolation to large radii does rely upon the faithful 
representation of an object by a shapelet expansion, and the re- 
moval of its noise via series truncation. Such truncation restricts 
the basis functions' completeness, and a weight function (con- 
structed from a combination of the allowed basis functions) akin 
to a "prior probability" is subtly implicit inside our fitting proce- 
dure. However, a fitting method like ours can beat a direct, pixel- 
by-pixel measurement. Our fit is able to include flux from the ex- 
tended wings of an object, by integrating it over a large area, even 
when the signal lies beneath the noise level in any individual pixel. 
The wings of galaxies in figure 7 are indeed well-captured by the 
shapelet models. 



Thanks to Mark Coffey for help deriving this expression. 



6.3 Size 

Measures for the size and ellipticity of an object can be derived 
from its unweighted quadrupole moments. 



Fij = / / XiXjf{-x.) d^x . 

The rms radius 71 of an object is given by 



(52) 



//Jx|V(x) d'x 



(53) 

^(n + l)/„o. (54) 



Integrals (46) to (49) can also be used to calculate Petrosian radii 
that enclose a specified fraction of the total flux within a circular 
aperture. 



6.4 Ellipticity 

The unweighted ellipticity of an object can also be calculated from 
its quadrupole moments. 



Fu - F22 + 2iFi2 _ (167r)5/3^ 



-Fii + F2: 



FR^ 



^ [n(n + 2)] i/„2, (55) 



where the complex ellipticity notation of Blandford et al. (1991), 
with e = |e| cos 20 -|- i|£| sin 20, arises here naturally. 



7 GALAXY MORPHOLOGY CLASSIFICATION 

The shapelet decomposition of an object captures its entire struc- 
ture, and useful information is frequently found in coefficients of 
higher order than those considered above. In particular, galaxy mor- 
phologies are well known to provide an indication of their physical 
properties, local environment and formation history. The classical 
"Hubble sequence" of morphological types has been recently im- 
proved by several shape estimators that attempt to classify galaxies 
in a more quantitative manner, which correlates directly with the 
physical properties of interest {e.g. Simard 1998; Bershady, Jan- 
gren & Conselice 2000; van den Bergh 2002). 

It is possible to manufacure such morphology diagnostics 
from weighted combinations of shapelet coefficients. Introducing 
shapelets to this field allows a measurement to take advantage of 
our robust treatment of noise, pixellisation and PSF deconvolution. 
The shapelet expressions for existing shape measures are frequently 
elegant; and the natural symmetries in shapelets also suggest new 
diagnostics. 
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7.1 General scale-invariant quantities 

One approach is to consider general shape estimators Q, formed 
from a linear combination of shapelet coefficients, as an extension 
of the previous section 



Q — ^ ^ ^n^mfrijrn ? 



(56) 



where arbitrary weights, and the exponent s sets the di- 

mension of the estimator. These are also linear in a galaxy's surface 
brightness. We initially restrict ourselves to using those combina- 
tions which are independent of /3 to at least first order. This ensures 
that the choice of the scale factor does not affect the final result, 
and is also equivalent to invariance under object dilations (40). We 
can then impose further constraints that the estimator must be inde- 
pendent of or linearly dependent upon the various other operations 
described in §5. Setting ^ = and using the result that 



dp " 2^ { V(" + "^ + 2)(«-f» + 2) /n+2,m 

- ^{n + m){n - m) /„-2,m j , (57) 
it is easy to show that we require 



2s 



\/(n -h m)(n — m) 



W„-2,m 



{n + m — 2){n — m — 2) 
(n + m)(n — m) 



(58) 



Notice that all quantities so formed mix coefficients with only one 
value of I m I . This can be chosen to give Q the desired properties 
under rotation. Any term on the right-hand side should be ignored 
if it refers to non-existent states with negative n. The normalisation 
of the first term in each series, Wm,m, is arbitrary: this can be set to 
ensure independence to changes of object flux. 

Setting (s, m) = (1, 0), (2, 1), (3, 0) and (3, 2) recovers the 
flux F, centroid Xc, rms square radius and ellipticity e, up to 
the normalisation factor of for the latter three quantities. This 
proves that these are indeed the only /3-invariant linear quantities 
with such dimensionality and rotational symmetries. Furthermore, 
since equations (50), (51), (54) and (55) describe for unweighted 
moments, they must in fact be independent of B to all orders. 

All of these basic shape estimators converge for any galaxy 
with a shapelet spectrum steeper than n~^. This includes both spi- 
ral galaxies with an exponential profile, and elliptical galaxies with 
a "de Vaucouleurs" profile, as long as rimax is kept sufficiently 
low to prevent the high-n coefficients from modelling background 
noise at large radii. The flux and centroid estimators converge most 
rapidly, so are least sensitive to the choice of nmax- The error on 
these series due to truncation can be calculated using any of a range 
of methods for generic Taylor series in e.g. Boas (1983). 



7.2 Concentration 

We can extend this sequence by raising s further. For example, set- 
ting s = 5 and m = gives the 2D unweighted kurtosis of the 
image, producing an estimate of the object's concentration. Unfor- 
tunately, such a high value of s yields a series of shapelet coeffi- 
cients that does not converge for galaxies with a de Vaucouleurs or 
exponential radial profile. 

We have also noticed that a ratio of the two existing shapelet 



scale sizes, R and 13, also works rather well as a concentration index 
(although this estimator is not independent of /?). Further work will 
need to be done to calibrate this estimator to the physical properties 
of galaxies. 

An alternative approach is to mimic pre-existing, and pre- 
calibrated, morphology diagnostics. Bershady, Jangren & Con- 
selice (2000) define a concentration index 



C s 5 X log 



r-za 



(59) 



where rgo and r2o are the radii of circular apertures containing 
80% and 20% of the object's total flux. This correlates well with 
a galaxy's Hubble type (Bershady, Jangren & Conselice 2000) and 
also its mass (Conselice, Gallagher, & Wyse 2002). Integrals (46) 
to (49) can be evaluated for various values of 7?, to find rso and 
r2o, and thus calculate this quantity for a shapelet model. 

7.3 Asymmetry 

Conselice, Bershady & Jangren (2000) define an asymmetry index 



A 



(60) 



Epixels /(^' y) 

where the superscript denotes an image rotated through 180°. A 
term dealing with the background noise and sky level has been 
omitted here, as these are automatically dealt with during the 
shapelet decomposition process in §4. Asymmetry correlates with a 
galaxy's star formation rate (ConseUce, Bershady & Jangren 2000), 
and high asymmetry values often indicate recent galaxy interac- 
tions or mergers. 

In a shapelet expansion, all of the information about galaxy 
asymmetry is contained in coefficients with odd m (and n). Using 
the orthonormality condition (9) and rotating via equation (38), we 
find the simple form 



A = 



V2/3 



(61) 



Estimators of asymmetry under rotations of 120° or 90° can 
also be formed from sums of shapelet coefficients with m not di- 
visible by 3 or 4 respectively. 

7.4 Chirality 

A quadratic combination of shapelet coefficients can be used to de- 
scribe the "chirality" or "handedness" of an object. One dimension- 
less estimator x\m\ can be formed for each value of |m|, to trace 
the relative rotation of those coefficients, with increasing n. This 
is roughly equivalent to tracing the rotation of a galaxy's isophotes 
with increasing radius. For example, the galaxy shown in figure 3 
has two prominent spiral arms that unwind in a clockwise sense, so 
it has a large, positive value of X2. 

We require that the chirality estimators should be invariant un- 
der global rotation of the object; invariant under changes of flux; 
invariant to first order under changes of (3; and to flip sign when the 
object is mirror-imaged. These conditions uniquely specify 



2 oo oo 



(62) 



X|77i| jp'2 ^ ^ ^ ^ ^n,n' fn,mfn' ,m 5 
n=m n'=n-\-2 

where 'Wm,m+2 = 1 and 

{n' + m)(n' - m)wn,n' = ^Wn,n'-2-\-V {n + m)(n - m) ,(63) 
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thus mixing all coefficients with the same value of m. 

This estimator has yet to be calibrated against physical quanti- 
ties. However, this approach ought to be able to automatically dis- 
tinguish between elliptical galaxies and spiral galaxies in way that 
mimics a visual classification, and could also be adapted as a func- 
tion of rimax to find bars in the cores of spiral galaxies. 



8 CONCLUSIONS 

We have extended the formalism of shapelets for image analysis 
from basis functions separable in Cartesian to polar coordinates. 
Cartesian shapelets are convenient for the initial object decompo- 
sition. In particular, we have shown that that they can anals^ically 
be integrated inside a square boundary, thus facilitating the pix- 
ellisation of the smooth basis functions. On the other hand, polar 
shapelets decompose an object into components with explicit rota- 
tional symmetries, and often have a more direct physical interpreta- 
tion. In addition, they yield more compact representations of typical 
galaxy images, since terms with low orders of rotational symmetry 
tend to dominate. 

We have quantitatively investigated the effects of the choice 
of the shapelet scale size parameter, j3. For most objects in astro- 
nomical images, one scale size is clearly optimal for high quality 
image reconstruction, data compression and the fast convergence of 
shape estimators. We have developed a practical algorithm to find 
this value of (3 for arbitrary objects in real images, plus optimimi 
values for the shapelet centre Xc and truncation order rimax. This 
algorithm can also take into account observational effects including 
noise, pixellisation and PSF deconvolution. 

We have then described a number of applications of po- 
lar shapelets. Shapelet models can be rotated, enlarged and 
sheared by simple analytic operations. Since the shapelet ba- 
sis functions are invariant under Fourier transform, analytic con- 
volutions and deconvolutions (e.g. from a PSF) are also easy 
to perform. Linear combinations of an object's polar shapelet 
coefficients produce elegant expressions for its flux (photome- 
try), position (astrometry) and size. We also showed how other 
combinations of shapelet coefficients can be used to distin- 
guish between morphological types of galaxies. A complete 
IDL software package to perform the image decomposition 
and shapelet image analysis is publicly available on the world 
wide web at http://www.astro.caltech.edu/~rjm/ 
shapelets/. 
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APPENDIX A: LAGUERRE POLYNOMIALS 

Different conventions have been used to define the Laguerre poly- 
nomials, especially before the 1960s. The p\ term is omitted from 
equation (6) in many older books, and caution must be observed 
with the resulting relations. Several useful recursion relations can 
be derived to simplify their calculation {e.g. Boas 1983, Ch.l2), 
which we gather here, using our convention, for future reference: 



Ll{x) 
LUx) 

Ll{x) 



I - x + q 

q — 1 — X 



(Al) 
(A2) 



2-1- ■ 



P 



V 



- ( 1 + ^) Ll_^{x) (A3) 

(A4) 

(A5) 



AL%{x) 
da; 

d.L%(x) 
da; 



= x-^{vLl{x)-{jp + q)Ll_^{x)) (A6) 

= -l^tXix) (Al) 

= Lp-i{x) (A8) 

= -Xii?. (A9) 



This paper has been produced using the Royal Astronomical Soci- 
ety/Blackwell Science style file. 



